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We report a new computational method based on the recursive Green's function technique for 
calculation of light propagation in photonic crystal structures. The advantage of this method in 
comparison to the conventional finite-difference time domain (FDTD) technique is that it computes 
Green's function of the photonic structure recursively by adding slice by slice on the basis of Dyson's 
equation. This eliminates the need for storage of the wave function in the whole structure, which 
obviously strongly relaxes the memory requirements and enhances the computational speed. The 
second advantage of this method is that it can easily account for the infinite extension of the 
structure both into air and into the space occupied by the photonic crystal by making use of 
the so-called "surface Green's functions". This eliminates the spurious solutions (often present in 
the conventional FDTD methods) related to e.g. waves reflected from the boundaries defining the 
computational domain. The developed method has been applied to study scattering and propagation 
of the electromagnetic waves in the photonic band-gap structures including cavities and waveguides. 
A particular attention has been paid to surface modes residing on a termination of a semi-infinite 
photonic crystal. We demonstrate that coupling of the surface states with incoming radiation may 
result in enhanced intensity of an electromagnetic field on the surface and very high Q factor of the 
surface state. This effect can be employed as an operational principle for surface-mode lasers and 
sensors. 

PACS numbers: 42.70.Qs, 41.20.Jb, 78.67.-n 



I. INTRODUCTION 



Optical microcavities and photonic crystals (PC) have 
received increased attention in recent years because of 
the promising prospects of applications in a future gen- 
eration of optical communication networks 1 ' 2 . Examples 
of successfully demonstrated devices include lasers, light 
emitting diodes, waveguides, add-drop filters, delay lines, 
and many other 3 . 

By far the most popular method for the theoreti- 
cal description of light propagation in these systems is 
the finite-difference time-domain method (FDTD) intro- 
duced by Yee 4 . The success of the FDTD method is due 
to its speed, flexibility and ease of computational storage 
requirements. The limitation of the FDTD technique is 
related to the fact that the computational domain is fi- 
nite. As the result, an injected pulse experiences spurious 
reflections from the domain boundaries, which leads to 
mixing between the incoming and reflected waves, fn 
order to overcome this limitation the so-called perfectly 
matched layer condition has been introduced 5 . However, 
even with this technique, a sizable part of the incoming 
flux can still be reflected back 6 . In many cases the sepa- 
ration of spurious pulses is essential for the interpretation 
of the results, and this separation can only be achieved 
by increase of a size of the computational domain 7 . This 
may lead to a prohibitive amount of computational work, 
because the stability of the FDTD algorithm requires a 
sufficiently small time step. 

The problem of the spurious reflections from the com- 
putational domain boundaries does not arise in the meth- 
ods based on the scattering matrix technique, where the 



incident and outgoing fields are related with the help of 
the scattering matrix 8-12 . Other approaches where the 
spurious reflections are avoided include e.g. a multiple 
multipole method 13 , and a Green's function method 14 
based on the analytical expression for the Green's func- 
tion for an empty space. The main objective of the 
present paper is to present a novel computational ap- 
proach based on the recursive Green's function technique 
that can account for an infinite extension of a photonic 
crystal. In this technique the Green's function of the pho- 
tonic structure is calculated recursively by adding slice 
by slice on the basis of Dyson's equation. In order to ac- 
count for the infinite extension of the structure both into 
air and into the space occupied by the photonic crystal 
we make use of the so-called "surface Green's functions" 
that propagate the electromagnetic fields into infinity. In 
this paper we present a method for calculation of the sur- 
face Green's functions both for the case of a semi-infinite 
homogeneous dielectrics, as well as for the case of a semi- 
infinite periodic structure (photonic crystal) . This makes 
it possible to apply the Green's function technique for in- 
vestigation of a variety of important structures including 
waveguides and cavities in infinite or semi-infinite pho- 
tonic crystals, as well as to study the effect of the sur- 
face states and the coupling of waveguide Bloch modes to 
the external radiation. Note that the recursive Green's 
function technique is widely used for quantum mechan- 
ical transport calculations 15-18 and is proven to be un- 
conditionally numerically stable for various discretization 
schemes. 

The article is organized as follows. In Section II we 
present a general formulation of the problem. A descrip- 
tion of the recursive Green's function technique is given 
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in Section III. This section also provides a recipe for the 
calculation of Bloch states in a periodic structure as well 
as the surface Green's function. Technical details of the 
calculations are given in Appendices A-C. Several exam- 
ples of the application of the developed method are given 
in Section IV. The conclusions are presented in Section 
V. 



II. GENERAL FORMULATION OF THE 
PROBLEM 



For the numerical solution, Eqs. (4)-(6) have to be dis- 
crctizcd, x, y — > mA, nA, where A is the grid step. Using 
the following discretization of the differential operators in 
Eqs. (5),(6) 19 , 

2 d df(x) 
d 2 

A -^^(x)f(x) — ► Cm+l/m+1 - 2£ TO / m + £m-l/m-l 

(7) 



We start with Maxwell's equations in two dimensions we arrive to the finite difference equation 



j-^V x {V x E(r)} = ^E(r) (1) 

V X X H(r> } " ? H<r) ' 

where r = xi + yj, V = J^i + ^j, e r (r) is the 
relative dielectric constant, and the electric and mag- 
netic field vectors E(r, t) = E(r) exp(— iuit), H(r, t) = 
H(r) exp(— iojt). If the dielectric constant e r (r) is inde- 
pendent on z, the Maxwell's equations decouple in two 
sets of equations for the TE modes (H z , E x ,E y ), 

d 1 d TT d 1 d TT u? TT n 
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and for the TM modes (E z , H x , H y ), 
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Let us rewrite the equations for H Zl E z (2), (3) in an 
operator form 2 



(4) 



where the Hermitian differential operator L and the func- 
tion / reads, 

, TT T d 1 d d 1 d 
TE modes: / = H z , L TE = — - — , 

Ox e r ox Oy e r Oy 

(5) 

i / d 2 d 2 \ i 



TM modes: / = y/e^E z , L T m 



(8) 



^m.m\n,n-\-lfm.n-\-\ ^m.m-.n.n—lfm.n — l I ) fm.ni 

\ c J 

where the coefficients v, u arc defined for the cases of TE 
and TM modes as follows 

TE modes: / TOj „ = H zm y : £ m ,„ = — - — , (9) 

V m ,n = £ m + i >n + Cm-|,n + £m,n+| + £ro,ra-i' 
Um,m+l;n,n = Cra+i.nJ u m,m— l;n,n = £ m _I j7 j, 
U m ,m;n,n+1 = Cn,n+ij u m,m;n,n — l = £m,n— ii 



TAI modes. f m .n \/ m,nE z m .n'i £m,n 



(10) 



^m,ra+l;nn £m,n£m-fl,n7 ^Tn,m—l;nn ~ Cm-l,nCm,ni 
U m .rrr,n,n+1 6fi,n+l^m.nj u m,m;n,n—l — £,m,n£,m,n — l- 

A convenient and common way to describe finite- 
difference equations on a numerical grid (lattice) is 
to introduce the corresponding tight-binding operator. 
For this purpose we first introduce creation and an- 
nihilation operators, a+ n , a m . n - Let the state |0) = 
|0, . . . , m ,n, ■ ■ ■ , 0) describe an empty lattice, and the 
state |0, . . . 0, l m ,n, 0, . . . , 0) describe an excitation at the 
site m,n. The operators a+ n , a m , n act on these states 
according to the rules 16 

a+jO) = |0,...0,l m ,„,0,...,0), (11) 
o+ n |0,...0,l m ,„,0,...,0) =0, 



and 



\dx 2 dy 2 J \fs^ 
(6) 



a m ,„|0> = 0, (12) 
a m ,„|0,...0,l m ,„,0,...,0) = |0). 
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and they obey the following commutational relations 

(13) 



[Q>m,n,Q'm,n\ — [ a rn.m a m'n 



0. 



Consider an operator equation 



t f {— ) />■ 
where the Hermitian operator 

m,n 



(14) 



(15) 



«m, m-\- l;n,n " , m,7i"'Ti+l)^ ^m+l,m;n,n6E< m _|-i n &m,n 



acts on the state 



!/> = £/» 



,n|0>- 



(16) 



Substituting the above expressions for C and |/) in Eq. 
(14), and using the commutation relations and the rules 
Eqs. (11)-(13), it is straightforward to demonstrate that 
the operator equation (14) is equivalent to the finite dif- 
ference equation (8). Note an apparent physical meaning 
of the last four terms in Eq. (15): terms 2 and 3 describe 
forward and backward hopping between two neighboring 
sites in the ^-direction, and terms 4 and 5 denote similar 
hopping in the y-direction. In the next section we outline 
the Green's function formalism for solution of Eq. (14). 



III. THE RECURSIVE GREEN'S FUNCTION 
TECHNIQUE 

A. Basics 
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FIG. 1: Schematic illustration of the system under study de- 
fined in a waveguide (supercell) of the width TV. An internal 
region of the structure occupies M slices. Two representative 
cases are shown, (a) external regions are semi-periodic pho- 
tonic crystals with the period M, (b) external regions repre- 
sent a semi-infinite periodic photonic crystal with the period 
Ai to the right and air to the left. Arrows indicate the direc- 
tions on the incoming (/), reflected (7?), and transmitted (T) 
waves. 



Let us hrst specify structures under investigation. We 
consider light propagation through a photonic structure 
defined in a waveguide (supercell) of the width N, where 
we assume the cyclic boundary condition (i.e. the row 
n = N + 1 coincides with the row n = 1). The photonic 
structure occupies a finite internal region consisting of M 
slices (1 < m < M). 

The external regions are semi-infinite waveguides (su- 
percells) extending into regions m < and m > At + 1 . 
The waveguides can represent air (or a material with a 
constant refractive index), or a periodic photonic crys- 
tal. Figure 1 shows two representative examples where 
(a) the semi-infinite waveguides represent a periodic pho- 
tonic crystal with the period M, and (b) a photonic 
structure is defined at the boundary between air and the 
semi-infinite photonic crystal. 

Let us first define the scattering states for the struc- 
tures under consideration. The translation invariance 



along the supercell dictates the Bloch form for the ath 
incoming state \ipa), 



N 

m<0 



n=l 



(17) 



where fc+ (k~) is the Bloch wave vector of the right- 
propagating (left-propagating) state a, and 4>m the 
corresponding Bloch transverse eigenfunction satisfying 
the Bloch condition 



J m,n Tm+M,n- 



(18) 



The transmitted and reflected states, |^) and \ip r a ), can 
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be written in a similar form, 



N 



K>= £ E** ett?(m_( " +1)) E<»<»i°>. 

/3 n=l 

(19) 

JV 

1^) ^ E E ^ ei ^ m E<n<n|0>, (20) 



m<0 /3 



n=l 



where tp a {rp a ) stands for the transmission (reflection) 
amplitude from the incoming Bloch state a to the trans- 
mitted (reflected) Bloch state (3. Note that in general case 
the wave vectors and the Bloch states 0° n can be dif- 
ferent in the left and right waveguides (see e.g. Fig. 1(a), 
when the photonic structure is defined at the boundary 
air/photonic crystal). The method of calculation of the 
Bloch states for an arbitrary periodic structure is de- 
scribed below in Section IIIC. 

We define Green's function of the operator £ in a 
standard way, 



(K/c) 2 - Z) G(u) = 1 , 



(21) 



where 1 is the unitary operator 20 . The knowledge of the 
Green's function allows one to calculate the transmission 
and reflection coefficients. Indeed, let us write down the 
solution of Eq. (14) as a sum of two terms, the incom- 
ing state IV' 1 ) and the system response \ip) representing 
whether the transmitted or reflected states, or W), 
|/) = + \ip). Substituting |/) into Eq. (14) and using 
the formal definition of the Green's function Eq. (21), 
the solution of Eq. (14) can be written in the form 



G 



(2 - (^A/c) 2 ) IV 1 ). (22) 



Calculating the matrix elements (M + l,n\tp) = 
(0|aM+i,nV') and {0,n\ip} = (0\ao, n '>P), of the rigth and 
left hand side of Eq. (22), we arrive to the N x N sys- 
tem of linear equations for the transmission and reflection 
amplitudes (see for details Appendix A), 

$ M+ iT = -G M+1 ' (C/ ,i$-M+i^ - rr^o) (23) 
$ R = -G°'°(Uoa^-m+iKi - rr^o) - $0 (24) 

where the matrix elements {T)p a = tp a , (R)p a = r f3a'i 
qm +i,o anc j qo,o are Q reen ' s function matrixes with 
the elements 



(G m '%, p = {0\a m , n Ga+ p |0>. 



(25) 



Ti = G^'g is the left "surface Green's function" corre- 
sponding only to part of the whole structure, namely, 
to the semi-infinite waveguide (supercell) that extends 
to the left, — oo < m < 0. The physical meaning of the 
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FIG. 2: Schematic illustration of the application of Dyson's 
equation for calculation of Green's function for a composed 
structure consisting of m+l slices (see text for details). The 

operators C m and I m+1 describe respectively the structure 
composed of m slices, and the (m + l)-th slice. The opera- 

tor Cm+i — C m + l m+ i + V corresponds to the composed 
structure ol m + l slices, where V is the perturbation opera- 
tor describing the hopping between the mth and (m + l)-th 
slices. 



surface Green's function V is that it propagates the elec- 
tromagnetic fields from the boundary slice of the semi- 
infinite waveguide (supercell) into infinity. A method for 
calculation of the surface Green's functions both for the 
case of a semi-infinite homogeneous dielectrics, as well as 
for the case of a semi-infinite photonic crystal in a waveg- 
uide geometry is described below in Section HID. The 
matrixes Ki and <I> m are given by the right-propagating 
Bloch eigenvectors fc+ and the corresponding eigenstates 
</>^ „ in the waveguides, 

{Ki) aj3 = exp(ik+) 5 a p; ($ m )„ Q = CU (26) 
and the diagonal "hopping matrix" Uq,i is defined as 



{Uo,l)n,n' — U0,l;n,n'8 n ,n'- 



(27) 



(Note that the matrix Ki in Eqs. (23), (24) refers to the 
right-propagating states in the left waveguide). In the 
following sections we describe the recursive Green's func- 
tion technique based on the successive use of the Dyson's 
equation, introduce the method for the calculation of 
Bloch states in a periodic structure, and outline the way 
to calculate the surface Green's function T. 



B. Recursive technique based on Dyson's equations 

In order to calculate Green's function of the internal 
structure (i.e. for the slices 1 < m < M) we utilize the 
recursive technique based on Dyson's equation, see Fig. 
2. 

In order to illustrate this technique let us consider a 
structure consisting of m slices. The operator C m de- 



■5 



scribing this structure can be written down in the form 
^rn = ^v r a+a r - ^u r!r+A a+a r+A , (28) 

r r,A 

where r = m',n' (1 < m! < to; 1 < n' < N), and the 
summation over A in the second term is performed over 
all available nearest neighbors. Suppose we know Green's 

function G„ of the operator £ m , as well as Green's func- 

tion g„ l+l of the the operator I m+1 corresponding to a 
single (m + l)-th slice, 

n 

^m+l,m+l;n,n+l^' m ^i. n O'7n+l,n+l 
u m+l .rn+l;n+l ,ri a m +i _ n +i&m+l .n) • 

(The method of calculation of Green's function for a sin- 
gle slice is outlined in Appendix C). Our aim is to calcu- 
late Green's function of the composed structure, G m +i, 
consisting of to + 1 slices. The operator corresponding to 
this structure can be written down in the form 

C m+1 = C° m +t m+1 +V, (30) 

where the operators C m and I m+1 are given by the ex- 
pressions Eqs. (30), (29), and V = V m ,m+i m+l.m 
is the perturbation operator describing the hopping be- 
tween the TOth and (to + l)-th slices, 

V = V m +l,m + V m,m+l = (31) 
n 

The Green's function of the composed structure, G m +i, 
can be calculated on the basis of Dyson's equation 20 

G m +i = G° + G°VG m+ i, (32) 
Gm+i = G° + G m +iVG°, 

where G° is the 'unperturbed' Green's function corre- 
ct) 

sponding to the operators C m or I m+1 - For the sake of 
completeness, a brief derivation of Dyson's equation is 
given in Appendix B. Thus, starting from Green's func- 
tion for the first slice c/° and adding recursively slice by 
slice we are in the position to calculate Green's function 
of the internal structure consisting of M slices. Explicit 
expressions following from Eqs. (32) and used for the 



recursive calculations are given below, 

n m+l,m+l _ ,j jj (r<0 \m,m TT \ — 1 „0 

^rn+l — y 1 ~ 9 m +l U m+l,m{tr m ) u m,m+l) 9 m+1 

(33) 

^ro+1,1 ^,m+l,m+lT7 \m,l 

/-yl,l (f-t0 i (fiO \\,mjj 

u m+l ~ V-'m) ' V-'m) l>m,m+l*J m +l ! 

^yl,m+l /^Q \l.mjj s-im+l,m+l 

^m+1 — V-'m) ( - / m,m+l l - r m +i > 

where the upper indexes define the matrix elements of 
the Green's function G m ' m ' = (0|a ro ,„Ga+,„,|0). This 
recursive technique is proven to be unconditially numer- 
ically stable 15-17 . The performance of the method is de- 
termined by the size of the system of linear equations 
(33) which we solve when we add each consecutive slice. 
This system is solved M times, where M is the number of 
slices of the internal structure (in the x-dircction) . The 
size of Eqs. (33) is N x N, where N is a number of dis- 
cretization points in the y-direction. Typical dimensions 
of the equations used for computations of the structures 
reported in Section 4 are ~ 200 x 200. 

In order to calculate the Green's function of the whole 
system, we have to connect the internal structure with 
the left and right semi-infinite waveguides. Starting with 
the left waveguide, we write 

£ int+left = £ int + £ left + V , (34) 

where the operators C int+left ■, £ int an d £ left describe 
respectively the system representing the internal struc- 
ture + the left waveguide, the internal structure, and the 
left waveguide. The perturbation operator V describes 
the hopping between the left waveguide and the internal 
structure. Applying then the Dyson equation in a similar 
way as we described above, 

Gint+left =G +G VGi n t+left, (35) 

we are in the position to find the Green's function 
Gint+left of the system representing the internal struc- 
ture + the left waveguide. G° in Eq. (35) in an 'un- 
perturbed' Green's function corresponding to the inter- 
nal structure and the semi-infinite waveguide (the "sur- 
face Green's function" T). Having calculated the Green's 
function Gint+left on the basis of Eq. (35), we proceed 
in a similar way by adding the right waveguide and cal- 
culating with the help of the Dyson's equation the total 
Green's function G of the whole system. 

C. Bloch states of the periodic structure 

In this section we describe the method for calculation 
of the Bloch states in periodic waveguides (supercells) 
using the Green's function technique. Similar method 
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pact form 



Ti 



4>M + 1 



= T 2 



V'o 



where 



(39) 



Ti = 



cell 
•^cell 



u- 



1,0 



To = 







^cell ^1,0 



with / being the unitary matrix. The wave function of 
the periodic structure has Bloch form, 



4>M+r, 



e ik * M Ii> m . 



(40) 



Combining Eqs. (39) and (40), we arrive to the 
eigenequation for Bloch wave vectors and Bloch states, 



FIG. 3: Schematic illustration of the calculation of Bloch 
states in an infinite periodic structure (see text for details). 
The operator £ C eii describes a unit cell under consideration, 
1 < m < A4, and £ ou t describes the rest of the structure. 
The hopping between the cell and the rest of the structure is 
described by the operator V . 



T7 X T 2 



ipo 



Ak x Ai 



ipo 



(41) 



determining the set of Bloch eigenvectors fc" and eigen- 
functions ip a , 1 < a < N. 

To improve numerical stability of Eq. (41), it may be 
rewritten in the form 11 : 



was used for calculation of Bloch states in quantum- 
mechanical structures 18 . 

Consider a unit cell of a periodic waveguide occupying 
M slices, 1 < m < M, see Fig. 3. 

Rewrite the operator corresponding to the whole struc- 
ture in the form 



V, 



(36) 



where the operators C ce n and C ou t describe respectively 
the cell under consideration (1 < m < JA), and the out- 
side region including all other slices — oo < m < and 
A4+1 < m < oo, and V is the hopping operator between 
the cell and slices m — and m = M + 1. Write the 
total wave function = J2 m ,n Vv«c 

\ip) = l^ccii) + |Vw) 



, n |0) in the form 
(37) 



where IV'ccii) and |"0out) are respectively wave functions 
in the cell and in the outside region. Substituting Eqs. 
(36), (37) into Eq. (14), we obtain |^ ceU ) = G ceU V |VWt), 
where G CG ii is the Green's function of the operator £ cc n. 
Calculating the matrix elements (l,n\ip) and (Ad^lip), 
this equation can be written in the matrix form, 



4>i 



'coll 
/-iM.Mrr I 

^ccll u i,oVm + i, 



(38a) 
(38b) 



where the vector column xp m = (ipm,i ■ ■ ■ VVn.Af) , and 
where we used Um,m+i — ^o,i (because of the period- 
icity) and C/o,i = Ui,o (according to the definition of U, 
Eq.(27)). It is convenient to rewrite Eq. (38a) in a com- 



(Ti+Ta)- 1 ^ 



ik x M 



1)" 



(42) 



This technique allows one to avoid overflows and un- 
derflows in the eigensolver routine when eigenvalues with 
| e *fc,A<| > i and | e *fc x Ai| ^ i are calculated. 

In order to separate the left- and right-propagating 
states we compute the Poynting vector integrated over 
transverse direction, whose sign determines the direction 
of propagation. Bloch state propagating in a waveguide 
(supercell) defined in a photonic crystal is illustrated be- 
low in Fig. 5(c). 

Poynting vector can be expressed as follows 2 



S a (y) = \m a {y) x H*(y)]. 



(43) 



Note that for the case of the waveguide defined in air, 
M = 1, and Green's functions G ce ii in Eq. (39) is sim- 
ply given by Green's function of a single slice g° (see 
Appendix C for details of calculation of g°). 



D. The surface Green's function T. 

Consider a semi-infinite Bloch waveguide (supercell) of 
the periodicity M. extending in the region — M < m < oo 
as depicted in Fig. 4 

Suppose that an excitation |s) is applied to its first slice 
m = —M. Introducing the Green function G wg corre- 
sponding to the operator C wg describing the waveguide, 
one can write down the response to the excitation \s) in 
the form 



IV-) = G wg |s), 



(44) 
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FIG. 4: A schematic diagram illustrating calculation of the 
surface Green's function T of a periodic structure (see text 
for details). 
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where is the wave function that has to satisfy Bloch 
conditions (40). Applying Dyson's equation between the 
slices and 1 we obtain 



G 



1,-M 
wg 



Ipt/i nG 



0,-M 
wg ) 



(45) 



where T r = Gg^ is the right surface Green's function. 
(Note that because the waveguide is infinitely long and 
periodic, Q% = G^ 1 ^ 1 = Glf+^ M +* = ... etc.). 
Taking the matrix elements (1, n\if>) of Eq. (44) and mak- 
ing use of Eq. (45), we obtain for an each Bloch state 
a, ip" = TrUi^tpo- The latter equation can be used for 
determination of T r , 



(46) 



where and are the square matrixes composed of 
matrix-columns tpf and ijiQ, Eq. (41). If the waveguide 
is open to the left, its surface Green's function is the 
same as the surface Green's function of the correspond- 
ing waveguide open to the right, Ti — T r . Note that 
for the case of the waveguide defined in air the surface 
Green's function (46) simplifies to T r U\$ = K, where K 
is defined according to Eq. (26). 



IV. APPLICATIONS OF THE METHOD 

To reveal the power of the method we study three 
model systems defined in 2D square-lattice photonic crys- 
tal. First, we calculate a transmission coefficient and 
quality factor (Q factor) of several representative types 
of microcavities in infinite PCs. Then we focus on semi- 
infinite crystals where we investigate the effect of surface 
states, and, finally, we consider a semi-infinite PC with 
a waveguide opening to the surface. For the bulk crys- 
tal we choose a structure composed of cylindrical rods 
with the permittivity e r — 8.9 and the diameter of a rod 
d = 0.4a in a vacuum background, where a is the size of 



FIG. 5: (color online) (a) Band diagram for the right- 
propagating TM-mode of an infinite 2D photonic crystal 
(e r = 8.9, d = 0.4a) in FX-direction. PC has a fundamen- 
tal bandgap in the frequency range 0.28 < uia/2irC < 0.44 
(filled with gray in the figure). Green line in the fundamental 
bandgap corresponds a guided mode in a waveguide created 
by removing a central row of rods from the PC as shown in 
the inset, (b) Additional bands (encircled with red) origi- 
nated from the finite size effect. The waveguide (supercell) 
contains three unit cells in the transverse direction as illus- 
trated in the inset, (c) Bloch state propagating in the PC 
waveguide at coa/2nc = 0.38 



the unit cell. Each unit cell is discretized into 25 points 
in both x and y directions. 

Most of photonic crystal devices operate in a bandgap. 
The structure at hand has a complete bandgap for TM- 
modes in the frequency range 0.32 < Loa/2nc < 0.44 1 , 
and does not have a complete bandgap for the TE- 
polarization. Because of this, we will hereafter consider 
the TM-modes only. 

The developed method allows one to treat structures 
unlimited in x-direction, whereas in y-direction the struc- 
ture of interest is confined within a supercell with im- 
posed cyclic boundary conditions. This leads to the fi- 
nite size effects in a photonic band structure. If the su- 
percell consists more than one elementary cell, additional 
bands appear along with the bands for infinite PC (Fig. 5 
(a,b)), as the result of the imposed boundary conditions 
in the transverse direction. 

A similar finite size effect emerges when air waveg- 
uides (supercells) are attached to the system of inter- 
est. Even though we send a wave from an open space, 
we use a finite number of propagating modes. Solu- 
tion of the eigenvalue problem (4) for the air super- 
cell give s a discrete set of r ight-propagating cigenstates 
k™ 1 = \Joj 2 I c 2 — (2irm/w) 2 where w is the width of the 
supercell, and m is integer such that max |m| < low/2-kc. 
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FIG. 6: Dispersion relation for the air supercell of the width 
of 9a. Effective angles of incidence are determined by the 
angular wavenumber m. Inset shows the effective angles of 
incidence for m = —2, — 1, 0, 1, 2. 



Thus, a wave incident from air effectively propagates only 
at certain incidence angles, determined by the ratio of the 
longitudinal and transverse wave vectors tan a = k™ / k™ , 
as illustrated in Fig. 6. Note that this finite size ef- 
fect (caused by the cyclic boundary conditions in the y- 
direction) might in some cases represent a drawback of 
the method. 



A. Microcavity 

In this section we consider a microcavity defined in 
a waveguide in an infinite PC. The waveguide is cre- 
ated by removing a single central row of cylinders, such 
that in the energy range corresponding the fundamental 
bandgap only one waveguide mode can propagate. Band 
diagram of the waveguide mode is shown in Fig. 5 (a). 

Three different cavities are introduced in order to show 
the effect of geometry and demonstrate the importance 
of proper design of a cavity. The first cavity is defined by 
two rods placed on the lattice sites, see insets in Fig. 7. 
In the second structure the diameter of the rods is dou- 
bled, and for the third cavity we place two rods from each 
side of the cavity to achieve better confinement. A de- 
pendence of the transmission coefficient on the incoming 
wave frequency is depicted in Fig. 7(a). We would like 
to stress that in the calculation of the transmission co- 
efficient, the incoming, transmitted and reflected states 
are the Bloch states of a waveguide (shown in Fig. 5(c)), 
such that all spurious reflections from PC interfaces or 
computational domain boundaries are avoided. 

The fundamental parameter of cavity resonances 
is their Q-factor defined as Q=2ttui* (stored en- 
ergy) /(energy lost per cycle), which can be rewritten in 
the following form: 



4 / S in dy 




FIG. 7: (color online) (a) Transmission coefficient of three 
cavity structures versus frequency, (b) Intensity of the E z - 
component of the electromagnetic field in the double- wall cav- 
ity at the resonance (ua/2-KC — 0.3952). 



where tt T M = J[ee \E z \ 2 + fi (\H x \ 2 + \H y \ 2 )]dxdy and 
Ote = / \po\H z | 2 + se (\E x \ 2 + \ E y \ 2 )]dxdy characterizes 
the energy stored in the system respectively for TM and 
TE polarizations and the integral over Si n is the incom- 
ing energy flux. Eq. (IV A) can be also expressed as a 
well-known relation Q = u>/ Auj where u> is the resonant 
frequency and Auj is the width of the resonant peak at 
half- maximum. 

The resonance peak for the single-wall cavity is cen- 
tered at u>a/2irc = 0.3952 and has Q factor 35.5. As 
expected, the highest Q factor (327.7) is achieved for 
the case of double-rod walls. Resonance peak in the 
case of larger rods is shifted to the higher energy val- 
ues (uja/2TTc — 0.4281) because of the decrease of the 
effective size of the cavity. The lower Q factor in this 
case (25.07) is because the larger rods disrupt destruc- 
tive interference in a bandgap of the PC. 

Note that the width of the supercell used in the com- 
putations has to be large enough to ensure that the in- 
tensity of the field decays to zero at the domain bound- 
aries. At the same time, it is desirable to have the size 
of the computational domain as small as possible. For 
the present computations, keeping this trade-off in mind, 
we have chosen a supercell consisting of 7 unit cells in 
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the y-direction. This choice seems to be sufficient, as the 
field intensity decreases by 5 orders of magnitude within 
the length of two lattice constants from the waveguide 
towards the supercell boundaries. 

Finally, to confirm our results and to verify the devel- 
oped method, we performed calculations for the cavities 
and waveguides in PC studied by Li et al. u and found a 
full agreement with their results. 



B. Surface states 

In the previous section we considered wave propagation 
in an infinite photonic crystal. Another aspect of inter- 
est is the effect of the surface in semi-infinite photonic 
crystals that can accommodate a localized state (surface 
mode) decaying both into air and into a space occupied 
by the photonic crystal 1,21 . In the present section we 
study the coupling between an incident radiation and the 
surface states. Note that a surface mode residing on the 
surface of an infinite (in the y-direction) photonic crystal 
represents a truly bound state with the infinite lifetime. 
However, because of the used cyclic boundary conditions, 
our system is effectively confined in the transverse direc- 
tion. As the result, the translation symmetry is broken, 
and the surface mode turns into a resonant state with a 
finite lifetime. Using the developed method, we calculate 
the Q factor of the surface modes. Our findings indicate 
that the surface modes, thanks to their high Q factors, 
can be used for lasing and sensing applications. 

We study two semi-infinite photonic crystal structures 
that support localized surface modes. In the first case 
a surface row of cylinders is composed of half-truncated 
rods 1 (structure 1), and in the second case the cylindri- 
cal and half-truncated rods in the surface row are inter- 
changed as shown in Fig. 8 (structure 2). In order to 
calculate the Q factor of the structures at hand, we illu- 
minate the semi-infinite photonic crystal by an incidence 
wave (that excites the surface modes) and compute the 
intensity of the field distribution. Note that the calcu- 
lated field distribution includes the contributions from 
both the surface mode exited by the incident light, as 
well as the incident and reflected waves. This leads to a 
nearly constant off-resonance background in the depen- 
dence Q = Q(oj) that is caused by the contribution of the 
incident and reflected waves in the total field intensity in 
Eq. (IV A) . To remove this background we calculate the 
Q factor of a structure without surface states. We choose 
this structure as a semi-infinite photonic crystal with all 
identical cylindrical rods, which is known not to support 
surface modes 1 . Then the obtained value is subtracted 
from the calculated value of the Q factor of the system 
under study. Note that in the calculation of the Q factor, 
the surface integration in Eq. (IV A) is performed over 
the area depicted in Fig. 8. 

Figure 9 shows a Q factor of structures 1 and 2 as a 
function of the frequency of the illuminating light. For 
both structures the Q factor reaches ~ 10 4 . Figures 8 



(a),(c) show i? z -field distribution for structures 1 and 2 
at the resonance. For a comparison, a field distribution 
for a structure that does not support a surface mode (a 
semi-infinite photonic crystal with all identical cylindri- 
cal rods) is shown in Figure 8 (e). In the latter case the 
field intensity rapidly decays into the bulk of the pho- 
tonic crystal, whereas for the structures supporting the 
surface modes, the intensity is strongly localized at the 
boundary row of rods. It is also worth to mention that 
for the latter case the intensity of the field in the surface 
mode exceeds the incoming light intensity by 4 orders of 
magnitude, such that the light intensity in the air region 
is not visible in the figures (compare 8 (a),(c) with (e)). 

One can easily estimate the position of the resonant 
frequency for the surface modes. Indeed, the outermost 
row of the cylinders (where the surface state resides) can 
be considered as a resonator with the characteristic reso- 
nant wavelengths following from the cyclic boundary con- 
ditions and given by X a = 2ir/k a , where 



a is the mode number and w is the width of the waveg- 
uide. The surface state for structure 1 exists only in a 
limited frequency interval, 0.33 < uja/2nc < 0.37 (the 
dispersion relation of the surface mode of this structure 
is given in Ref. 1). It follows from this dispersion relation 
that all the modes given by Eq. (47), except a = 4, are 
situated outside this interval, whereas the mode a = 4 
corresponds to the frequency wa/27rc = 0.365. This es- 
timated frequency agrees very well with the actual cal- 
culated resonant frequency ua/2nc w 0.359, see Fig. 9. 
Note that the calculated field distribution, Fig. 8 (a), is 
fully consistent with the expected pattern for 4th mode in 
the system of n = 9 cylinders. (This field distribution is 
determined by the overlap of the eigenstate correspond- 
ing to the eigenfrequency (47) with the actual positions 
of the cylinders in the outermost row) . Note that we per- 
formed calculations for different numbers of cylinders in 
the transverse directions [n = 5, ...,11), and we always 
find an excellent agreement with the predicted value of 
the resonant frequency uj a . 

Figures 8 (b),(d) show Poynting vector distribution 
for both structures at the resonance. For the struc- 
ture 1 the Poynting vector is "curling" along the bound- 
ary, showing a low speed of the surface state. In con- 
trast, for the structure 2, the Poynting vector exhibits a 
rapid flow of energy along the boundary. Another differ- 
ence between these structures is a very broad and rather 
strong "background" peak in the structure 2 in the region 
0.34 < uja/2nc < 0.35 (with Q factor up to - 100). The 
presence of such the peak indicates that the correspond- 
ing surface state can be rather robust to various kinds 
of imperfections that are always present in real struc- 
tures and which are known to broaden the resonances 
and lead to decrease of the Q factor 12 . These two ex- 
amples of photonic crystals illustrate, that with proper 
structure design one can engineer and tailor properties of 



10 



9 o o o 
#000 

#000 

#000 

D O O O 

I o o o 

• 000 
#000 

• 000 



(a) 



D O 


O 




• O 


O 


O 




O 


O 


5 


O 


O 


• O 


O 


O 


S . 


O 


O 


<2> 


O 


O 





O 


O 


D O 


O 


O 





• o o c 

• 00c 

• 00c 

• 00c 

• 00c 

• 00c 

• 00c 

• 00c 

• 00c 



0016 
0015 
1 0014 
0013 
0012 
O01i 
0.01 

tow 
ooe 

0007 

00D6 

0005 

0004 
O0O3 
OOOZ 
I 0001 



(e) 



FIG. 8: (color online) E z field and Poynting vector distributions for the structure 1 ((a),(b)) and for the structure 2 ((c), (d)) 
at the resonant frequencies (marked by arrows in Fig. 9). (e) E z field distributions for the structure that does not support 
surface modes (a semi-infinite photonic crystal with all identical cylindrical rods). In all cases the structures are illuminated 
by the incoming wave with the incidence angle a = arctan k y / k x — 34.7°. 
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FIG. 9: Dependencies Q = Q(ui) for structures 1 and 2 (solid 
and dashed lines respectively). Arrows indicate the reso- 
nances for which the field intensities and Poynting vectors 
are visualized in Fig. 8 



the surface states into the required needs. 

High values of the Q factors of the surface modes re- 
siding at the interface of the photonic crystal structures 
indicate that these systems can be used for lasing and 
sensing applications. The lasing effect has been demon- 
strated for different photonic crystal structures includ- 
ing band-gap defect mode lasers 22 , distributed feedback 
lasers 23 , and bandedge lasers 24 . Utilization of the high- 
Q factor of the surface modes represent a novel way to 



sustain lasing emission. To achieve lasing effect care- 
ful design of the surface and surface mode engineering 
should be performed and the developed method seems to 
be a suitable tool for this purpose. A detailed study of 
the surface modes for various surface terminations, their 
Q-values, and dispersion relations will be reported else- 
where. 



C. Waveguide coupled to the open space 

The last example of application of the method pre- 
sented here is a semi-infinite photonic crystal with a 
waveguide coupled to the surface, see Fig. 10. It has been 
recently demonstrated that a surface of a photonic crystal 
can serve as a kind of antenna to beam the light emit- 
ted from the waveguide in a single direction 25,26 . These 
findings outline the importance of investigation of the 
surface modes in the photonic band-gap structures that 
can eventually open up the possibilities to integrate such 
the devices with conventional fiber optic devices. 

In the present section we consider two different crys- 
tal terminations to illustrate the effect of the surface on 
propagation of the light emitted from the waveguide. In 
the first case the surface is composed of cylinders with 
parameters identical to those in the bulk of the crystal, 
and in the second case the surface cylinders are two times 
smaller than the cylinders in the bulk. 

The Bloch state propagating in a waveguide in the 
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(a) (b) 



FIG. 10: (color online) E z field distributions at the surface 
of a truncated photonic crystal with a waveguide, (a) The 
surface is composed of cylinders with parameters identical to 
those in the bulk of the crystal, and (b) the surface cylinders 
are two times smaller than the cylinders in the bulk. 



be employed as an operational principle for surface-mode 
lasers and sensors. 
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APPENDIX A: CALCULATION OF THE 
TRANSMISSION COEFFICIENT 

In this appendix we provide a detailed derivation of 
Eqs. (23), (24). Consider first an infinite periodic struc- 
ture in a waveguide (supersell) geometry. The a-th Bloch 
state in the lattice can be written in the form 



photonic crystal couples with the states in air and the 
resulting field distributions is shown in Fig. 10. The first 
structure does not support the surface mode, and hence 
the light intensity distribution in the air region exhibits 
a typical diffraction pattern. However, for the case of the 
second structure the field distribution in the air region is 
drastically different. In this case the Bloch state in the 
waveguide couples with the surface state localized at the 
crystal termination, such that the whole surface acts as 
a source of radiation. 



V. CONCLUSIONS 

We have developed a method based on the recursive 
Green's function technique for the numerical study of 
photonic crystal structures. The method is proven to be 
an effective and numerically stable tool for design and 
simulation of both infinite photonic crystals and pho- 
tonic crystals with boundaries. In the present method 
the Green's function of the photonic structure is calcu- 
lated recursively by adding slice by slice on the basis 
of Dyson's equation. In order to account for the infi- 
nite extension of the structure both into air and into the 
space occupied by the photonic crystal we make use of 
the so-called "surface Green's functions" that propagate 
the electromagnetic fields into infinity. This eliminates 
the spurious solutions (often present in the conventional 
FDTD methods) related to e.g. waves reflected from the 
boundaries defining the computational domain. The de- 
veloped method has been applied to scattering and prop- 
agation of electromagnetic waves in photonic band-gap 
structures including cavities and waveguides. In partic- 
ularly, we have shown that coupling of the surface states 
with incoming radiation may result in enhanced inten- 
sity of the electromagnetic field on the termination of 
the photonic crystal and very high Q-factor of the sur- 
face modes localized at this termination. This effect can 



l^) = E e<ro ^,»<n|0), 



(Al) 



where summation is performed over all lattice sites and 
the function (j)^ nn satisfies the conditions (18). Substi- 
tuting Eq. (Al) into Eq. (14), we arrive to the finite 
difference equation valid for all sites m, n 



Vm-l,n 

(A2) 



-u 



Consider now the incoming state \ip l a ), Eq. (17). Sub- 
stituting Eq. (17) into Eq. (14) and using Eq. (A2) we 
obtain 



(A3) 



X] M l.O;™,n0o,„aJ n |O). 



Substituting this equation into Eq. (22), calculating the 
matrix elements (M + 1, n\tp) and (0, n\ip) , and using the 
relations 



Or — — (_r UiflL I, 

G o,o =Ti _ G0 .i u T 



(A4) 
(A5) 



that follow from Dyson's equation, we arrive to Eqs. 
(23), (24) determinig the transmission and reflection am- 
plitudes. 



APPENDIX B: DERIVATION OF THE DYSON'S 
EQUATION 

Let £° be the operator describing an unperturbed sys- 
tem and V be a perturbation. In our case the unper- 
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turbed system consists of several subsystems, e.g. m 
slices of the internal structure and (m + l)th slice, and 
the perturbation corresponds to the coupling (hopping) 
between them (see Fig. 2). The operator of the total 
(perturbed) system reads 



£° + V. 



(Bl) 



Let G° and G be the Green's functions of the unper- 
turbed and the total (perturbed) systems respectively. 
Starting with the definition of the Green's function (21), 
we obtain 

G- 1 = {uA/cf -£ = (luA/c) 2 - £° - V = (G )- 1 - V . 

(B2) 

Multiplying this expression from the left with G and from 
the right with G° we arrive to Dyson's equations 

G°G~ 1 G = G°(G°)~ 1 G - G°V G => 

G = G° + G°VG. (B3) 

Similarly one can also show that 

G = G° + GVG°. (B4) 



N 



(CI) 



n=l 



Using this operator in the definition of Green's function 
(21), and calculating the matrix elements {■■■) m ,m;n,ri = 
(0|a mj „... , |0), we arrive to the NxN system of linear 
equations for the matrix elements of the Green's function 
of a single slice g m , 



N / / a \ 2 

ujA 



n" = l \ V 7 / 

where the matrix element of the operator l m reads, 



m.m;n,n" J 9m,m:n" ,n' — &n,n' i 

(C2) 



— U m ,rn:n" -l,n" 5 n ,n" -1 — Um,m;n"+l,n"^n,n"+l- (C4) 



APPENDIX C: THE GREEN'S FUNCTION FOR 
A SINGLE SLICE 



The operator describing the m-th slice has the form 



Note that because of the cyclic boundary conditions 
in the /^-direction, the matrix elements u m rn: i N 

and 

u m ,m:N,i are distinct from zero and defined according to 

U m ,m;N,l — u m,m;0.1, and U m ,m;l,jV = Um,m;N+l,N ■ 
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